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ABSTRACT 

We consider the hydrodynamics of the outer core of a neutron star under conditions when both neutrons 
and protons are superfluid. Starting from the equation of motion for the phases of the wave functions of the 
condensates of neutron pairs and proton pairs we derive the generalization of the Euler equation for a one- 
component fluid. These equations are supplemented by the conditions for conservation of neutron number 
and proton number. Of particular interest is the effect of entrainment, the fact that the current of one nucleon 
species depends on the momenta per nucleon of both condensates. We find that the nonlinear terms in the 
Euler-like equation contain contributions that have not always been taken into account in previous applications 
of superfluid hydrodynamics. We apply the formalism to determine the frequency of oscillations about a state 
with stationary condensates and states with a spatially uniform counterflow of neutrons and protons. The 
velocities of the coupled sound-like modes of neutrons and protons are calculated from properties of uniform 
neutron star matter evaluated on the basis of chiral effective field theory. We also derive the condition for the 
two-stream instability to occur. 

Subject headings: hydrodynamics — stars: interiors — stars: neutron — pulsars: general 


1. INTRODUCTION 


The outer core of a neutron star consists of a uniform fluid 
of neutrons, protons and electrons, with possibly other minor¬ 
ity constituents. The hydrodynamics of the core of a neutron 
star is important for studies of a variety of phenomen a, amo ng 
them stellar oscillations (|Mendell|l99l||Lindblom & Mendel!] 


|1994||2000[|Andersson & Comer|200lji, co llective modes of 

matter (|Epstein||1988| |Bedaque & Reddy||2014|», as well as 
theories of spin-down and glitches in the rotation rate of neu¬ 
tron stars (For a review, see Haskell & Melatos~] ( |2015| )). From 
microscopic calculations, protons are expected to be super¬ 
conducting in the outer core, while the situation for neutrons 
is less clear because of the difficulty of calculating superfluid 
gaps at such densities with confidence. In this paper we shall 
consider the case when the protons are superconducting and 
the neutrons superfluid. 

The purpose of this paper is to derive the equations govern¬ 
ing the long-wavelength, low-frequency behavior of the sys¬ 
tem. We shall assume that thermal effects may be neglected: 


typical temperatures in neutron stars are of order 10 S K or 
10 keV, which is small compared with the Fermi energies of 
the components, which are of order 100 MeV. We shall fur¬ 
ther assume that the superfluid gaps are large compared with 
the thermal energy kgT, where kg is the Boltzmann constant 
and T the temperature. The basic variables in the approach 
we shall adopt are the density of neutrons, the density of pro¬ 
tons, and the phases of the condensate wave functions of pairs 
of neutrons and pairs of protons. This leads naturally to a 
description of the dynamics in terms of the gradients of the 
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phases, which correspond to the momentum per particle of the 
condensates. We shall show that this approach leads straight¬ 
forwardly to equations for the dynamics, including nonlinear 
terms, which agree with the work of Mendell ((1991). 

Of particular interest in this paper is the influence of en¬ 
trainment, the fact that there is a coupling between the cur¬ 
rents of the two components. To make the exposition as clear 
as possible, we shall derive the equations of motion by pedes¬ 
trian methods. We shall then show how they may be obtained 
from a Hamiltonian approach that exploits the fact that the 
phase of the condensate wave function of a component is the 
canonically_conjugate variable to the density of that compo¬ 
nent ( |Lifshitz & Pitaevskii||1980) . A particular focus of the 
work is to generalize the Euler equation for a one-component 
fluid to a two-component system, and we shall show that, in 
the Euler equations, there are contributions in the nonlinear 
terms in the Euler-like equations that have not always been 
considered in past applications, although they are implicit in 
the basic formalism (see, e,g, [Mendell| ( |199l) ). These arise 
because the quantity determining the degree of entrainment is 
a function of the densities of the two components. A prelimi¬ 
nary report of many of the results in this article was given by 


Kobyakov et al. (2015). 

This" article is arranged as follows. The basic formalism is 
described in Section 2, where we work in terms of the phases 
of the wave functions for the neutron and proton pair conden¬ 
sates, and the neutron and proton number densities. The equa¬ 
tions of motion for the phases are described by a Josephson 
equation and that of the nucleon densities by continuity equa¬ 
tions. Because of entrainment, the neutron number current 
depends not only on the gradient of the phase of the neutron 
condensate but also on the phase of the proton condensate. 
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and similarly for the proton current. In addition, entrainment 
affects the chemical potentials of nucleons. The specific form 
of the Euler-like equations for the momentum per particle of 
the condensates is derived in Section 3. Collective modes of 
oscillation about an initial situation where the condensates are 
stationary and the densities uniform are considered in Section 
4. There we also small deviations from a state with a uni¬ 
form counterflow of neutrons and protons. Applications to the 
outer core of a neutron star are described in Section 5, where 
we calculate collective mode velocities. Section 6 contains 
a general discussion of, among other things, the relationship 
between our work and some earlier work on superfluid hydro¬ 
dynamics. 


2. BASIC FORMALISM 

We shall consider long-wavelength, low-frequency phe¬ 
nomena, in which local charge neutrality is maintained and 
electrical currents are absent. This is a good approximation 
for frequencies small compared with the electron plasma fre¬ 
quency and for wavelengths long compared with the Debye 
screening length for electrons. Moreover, the hydrodynamic 
approximation implies that the frequency is smaller than 
the inverse of the electron relaxation time due to electron- 
electron collisions. We shall also neglect dissipation due to 
Landau damping of the electron motion, which will be treated 
elsewhere (|Kobyakov et al.|2016jl. Under these conditions, the 
system behaves as a two-component system, one component 
being the neutrons and the other the protons and electrons. 
Throughout we shall work in an inertial frame of reference, 
and therefore the centrifugal and Coriolis forces will not ap¬ 
pear explicitly. We denote the phase of the superfluid order 
parameter for neutrons by 2 (j>„ and that for protons by 2 (j) p . To 
first order in the gradients of the phases, one may write the 
number current density of neutrons as 


for protons. A separate continuity equation for electrons is 
not required, since the electron number density and current 
density are the same as those of the protons. 

We are interested in situations where spatial variations are 
slow. To determine how the phase of a state varies in time, we 
may therefore consider states in which the densities of neu¬ 
trons and protons are uniform, and the gradients of the phases 
are uniform. The equations of motion for the phases may be 
obtained by making use of the fact that in a state with energy 
£ the wave function varies in time as e | n terms of the 

ground states |(V„, N p ) of the system with N p protons and N n 
neutrons, the superfluid order parameter for neutrons is 

(Nn 2, N p | (r) y/ f/ |(r) \N n . Np) 

^ e -i [6[N„, N p )-S{N n - 2, N p )]t/h ^ 


where y/ na ( r) is the annihilation operator for a neutron of spin 
crj^l We remark that the energies of the states are also func¬ 
tions of the gradients of phases of the superfluid order param¬ 
eters. The quantity £(N n , N p ) — £(N„ — 2, N p ) is twice the 
neutron chemical potential including the contribution due to 
motion of the components, which we denote by JLi« ot Q From 
Eq. (|5]i we conclude that 


tot 

dt ’ 


(6) 


which is essentially Josephson’s equation (see, e.g., 
Varaquaux (|2015|l). In this article, we shall include in 


the calculations terms of second order in p a , and therefore 
we need the Hamiltonian to this order. In the Hamiltonian 
formalism, the current density is given by 


J a — 


8JP 

<5p« ’ 


(7) 


l np 


p P 


and that for protons as 


ip — 


l pp 


P p- 


L pn 


(1) 

( 2 ) 


where p a = 7zV0 a is the momentum per particle associated 
with the condensate and the response functions n n p .{a. p = 
n,p} generally depend on the density of neutrons, n n , and the 
density of protons, n p , but are independent of the gradients of 
the phases. Which mass is inserted in these equations is arbi¬ 
trary, but the choice of the nucleon mass m makes for simple 
expressions later in the analysis. To avoid inessential compli¬ 
cations, we shall neglect the difference between the neutron 
and proton masses. The quantity n a p/m is the long wave¬ 
length limit of the zero frequency neutron-current-density- 
proton-current-density response function and it is symmetri¬ 
cal in the indices a and /3. 

We shall assume that characteristic times for weak inter¬ 
action processes are long compared with the timescales of the 
motions, and therefore the numbers of neutrons and of protons 
are separately conserved. The continuity equation is therefore 


dn„ 

dt 

for neutrons and 

(9/7,- 


f ftnn . ftnp 

-P n + —P p 

v m m 


__L + V .(^ P; , + ^p„ 
dt \ m m 


= 0 


= 0 


(3) 

(4) 


where Ji? - f d?r H is the Hamiltonian, H being the Hamil¬ 
tonian density. 

It follows from Eqs. ([!} and ([ 2 J 1 that the “kinetic’]^] contri¬ 
bution to the Hamiltonian density i^] 


H km 


2m I " 2 m 




L np 


Pn * P p 


(ftnn T ftnp ) 2 

~^r^ Pn+ 


(jlpp T tlpn 

2m 


) 2 
Pp 


Hnp 

2m 


( 8 ) 

(Ph-Pp) 2 - 

(9) 


The final term in Eq. (|9]> represents the effects of entrainment 
of the motions of the two fluids. 

In the hydrodynamic description of a one-component fluid, 
the quantity <$> = h(j)/m is commonly referred to as the veloc¬ 
ity potential, since the fluid velocity is Vd>. However, we see 
from the considerations above that in multi-component sys¬ 
tems, the phases are more properly regarded as momentum 


3 For simplicity we consider the case of an S-wave superfluid, where pair¬ 
ing is in a spin singlet state. For superfluids with anisotropic gaps the pairing 
amplitude must be defined for particles with specified momenta. 

4 Since the phase is proportional to the difference of the energies of ground 
states whose neutron numbers differ by 2, it is a smooth function of N n and 
does not depend on whether N n is odd or even. 

5 We shall refer to as “kinetic” all contributions due to the motion of the 
components, including that due to entrainment. 

6 In the literature, the symbol v is used to denote an average velocity in 
some places and the momentum per unit mass of the condensate particles in 
others. To avoid confusion, we shall generally work with p, the momentum 
per particle in the condensate. 
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potentials, since the momentum per particle of species a in 
the condensate is TzV0 a . 

The remaining contribution to the Hamiltonian density is 
the energy density of the system in the absence of gradients 
of the phases, which we denote by E(n„,n p ). Thus the Hamil¬ 
tonian density is 


jkin 


H = E(n„,n p )+H k 

The equations of motion for the phases are therefore 

8J4? 


( 10 ) 


h -j— = — 
dt 


8n„ dn„ 


dH dH km 

= - 


and 


where 


dt 


f~tn — 


SJ4? _ 

8n p 

dE(n„,n p ) 

dn„ 


dH 

dn„ 


Pp 


dn„ 

d H kin 
dn n 


and p p = 


dE(n„,n p ) 

dn„ 


(13) 


are the neutron and proton chemical potentials when the 
phases of the condensates do not vary in spaceQ Since we 
consider matter that is electrically neutral, the quantity fj. p is 
the energy to add an electron and a proton but, for notational 
simplicity, we do not indicate this explicitly. 

Quite generally, the equations of continuity for neutrons 
and for protons have the form 


and 


The neutron density and the phase of the neutron conden¬ 
sate are conjugate variables, and these results also follow from 
the Hamilton equation, dn„/dt = with the expres¬ 

sion for the current given in Eq. (Ilbpj Similar results hold for 
the protons. In the Hamiltonian formalism, the “coordinates” 
and “momenta” are to be regarded as independent variables. 
Consequently, the derivatives on the right hand sides of Eqs. 
(JTTJ) and (JT2jl are to be evaluated at fixed p„ and p p . 

ine basic thermodynamic identity at zero temperature may 
thus be written as 


dn„ _ . 


¥ +v ' J " = 0 ' 

(14) 

dn„ _ . 

DT + V ^ = °’ 

(15) 


dE tot = jl^dtin + jl { p { dn p + j„ • dp n + j p • dp p , 


(16) 


where the energy density E lot and the chemical potentials ji P / 
all include kinetic contributions. 

The velocities of the components are defined by 


J a 

Va = —, 

n a 


(17) 


and therefore it follows from Eqs. ([TJ and ([2]) that 


1 


V«= -P|B- < 18 ) 

n a m p 

7 From the discussion after Eq. it follows that the deriva¬ 
tive dE/dn„ must be regarded as the limit for small integer V of 
[£(N n ,Np) — S(N n — 2v,N p )\/2vV, where V is the volume of the system. 
Similar results apply for the proton chemical potential. Consequently, odd- 
even effects due to pair breaking do not enter in the derivatives. 

8 Strictly speaking, the conjugate variables are ti<f) a and n a but we shall 
generally work in units in which ft is equal to unity. 


For a Galilean-invariant system, there are relationships^ be¬ 
tween the n a p (Mendell||l991 [Borumand, Joynt, & Kluzniak 
1996). Under a transformation to a frame moving with re¬ 
spect to the original frame by a velocity —v, the phases 0 a are 
increased by an amount m\ ■v/Ti. Consequently, the current 

)v- 


L np 


density of neutrons is increased by an amount (n„ 

However, from Galilean invariance, the change in the neutron 
current density is n„\ and therefore 



Wnn H - tlnp — fin- 

(19) 

(11) 

Similarly, by considering the proton current density one finds 


n pp T n np — n p . 

(20) 


Therefore Eq. (|9]l may be written as 


(12) 

rrkin n n 2 , n P 2 n "P ( n „ ^2 

H -2m P " + 2m Pp -2m (P "- Pp) - 

(21) 


3. EULER EQUATIONS 

The generalizations of Euler’s equation for a single com¬ 
ponent fluid to the two-fluid case are obtained by taking the 
gradient of Eqs. <E} and ([T2| and have the form 


dpn 

dt 


= -V jU„ 


Pi 


1 dn 


np 


2m 2m di 


(P»-Pp) : 


and 


dp P 

dt 


= - v AV 


P 2 p 


1 dn 


np 


2m 2m dt 


(P»-Pp ) 2 


( 22 ) 


(23) 


since p t/ — V 0„ ■ We may write the terms nonlinear in the p a 
in Eqs. (22 1 and (23 1 by using the vector identity V (a 2 /2) = 
a-V a + axV xa. Since in this article we shall consider 
only situations in which there are no singularities in the flow, 
we may put V X p n = 0 everywhere, and therefore 

dpn 1 — 1 8n np , , 

— + -p„ ■ Vp„ - - (p„ - p p ) • V(p„ - p p ) 

dt m m dn n 


— — VPn + — 

2m 


T^ dtlnp 


(P«-Pp ) 2 


(24) 


and 


1 dn 


dp P 1 „ 

dt m m dn 


"^ VMp+ 2 m 


,,Pt Vn-Vp)- V(p„-Pp) 

Pn-Vp) 2 - 


np 


(25) 


An interesting point is that the additional contribution to the 
nonlinear terms in the generalization of Euler’s equation is 
proportional to derivatives of n np , as feature already present 
in the work of Mendell (1991 Eqs. 14, 15, 29 and 30). 

Equations (|2?]l and <|25[) may be expressed in terms of 
the velocities or the components, but the resulting equations 
are lengthy because of the numerous places where density 
derivatives of n np appear. In the Appendix we show that the 
Euler-like equations in some earlier work do not agree with 
Eqs. ([22} and (23 1 . 


4. COLLECTIVE MODES 
4.1. Linear modes 
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We first consider the frequencies of modes corresponding 
to small deviations from the situation in which both superflu¬ 
ids are at rest (p„ = p p = 0). For a disturbance ~ c l!fc ' r m \ 
the perturbations of p a must be in the direction of k and 
Eqs. ( p~4| >, (p~5|>, (f24[ and (251 when linearized may be writ¬ 


ten in the matrix form 


—mv 

0 

ftnn 

Mfip 

0 

—mv 

Wpn 

n pp 

Enn 

Enp 

—V 

0 

E np 

E pp 

0 

— V 



(26) 


where v = <x>/k is the phase velocity of the wave. The mode 
frequencies are determined from the zeros of the determinant 
of the matrix, i.e.. 


v 4 - c 2 s v 2 + -4 det[£ ajS ] det[n ap ] = 0 , 


where 


det[E ap ]det[n a p] 


-'np ,l np 



C? = 


J PP' 1 PP 


(27) 


(28) 


(29) 


Equation (28 i is a generalization of the result of Bedaque & 
Reddy (|20T4f to allow for entrainment. In the absence of cou¬ 


pling between neutrons and the charged particles ( E np = 0 , 
= 0 , n, m = n„, and n pp = n p ), the mode velocities are given 


L np - 


by 


v 2 = 


for the neutrons and 


v° = 
p 


PnJ’-n 


P-p Ep p 


1/2 


1/2 


(30) 


(31) 


for the charged particles. 

One sees from Eq. ( |28] > that mode frequencies become 
imaginary if det[E a p] or det[n a p] become negative. The 
first condition corresponds to an instability to formation of 
a density wave with proton and neutron densities in phase 
for E np < 0 and out of phase for E„ p > 0. Generaliza¬ 
tions of this result to finite wavelengths have previously 


been employed to obtain estimates of the density at which 
the transition between uniform matter at higher densities 
and inhomogeneous matter in the crust occurs (|Baym,| 
|Bethe, & Pethick||1971[ |Hebeler et al.| 201 3| >. The con¬ 
dition det[n a p\ < 0 signals an instability to counter flow 


of the two components but, as we shall see in Section 5.3 


this is not expected to occur in the outer core of a neutron star. 


4.2. Two-stream instability 

We now consider small perturbations about a state in which 
the densities are uniform and the gradients of the phases are 
also uniform with values pj] and p° p . On linearizing Eqs. (14 1 , 
( fT5) , ( [34] ) and ( [35] ), one finds 


dn n n nn ~ n np 

- 5 - H-V • <5p„ + -4- V • Sp p 

at m m 

IjO p^ 

+--v«„„ + ^-v« 0 , 

m m 

dn p n pp n np 

- 37 - + — V • 5p p -|- L V • <5p„ 

at 111 m 

P° D° 

+ ^-Vn VD +^-Vn„ D = 0 , 


L PP ' 

m m 


(32) 


(33) 


<?c5p„ 


, — • V5p„ - — ^^( p 2 - P° p ) • V(5p„- <5p p ) 
dt m m dn n ‘ 


vk- 


1 dn, 


n P /_ 0 


2 111 dn n 


(Pii-Pp) = 0, (34) 


+ . V 5 p p - -^lT(p 0 _ p 0 ). V ( gpn _ gp ) 

dt m m dn r 1 


<9<5p 


v [Up - 


1 


L p 

dn, 


np / o 


2m dn u 


(P2-Pp) 2 )=0, (35) 


where 5p„ = p„ — pJJ and 5p p = p p p° On physical grounds 
one expects the most unstable mode to be one in which the 
wave number, and therefore also the velocity perturbations, 
are parallel to Ap = pjj — p°. For that case, Eqs. (32 1 , (331, 
([34| and ([35]) may be written in the matrix form 


m(0 

k 


F — 

t-'nn 


Enp 


+ P°n~^ 



n PP 

ftnp 

-*p d iz 

m(0 

k 

+ Pp~ A P^ 

ftpn 

^nn 

1 d2n " p (\n) 2 

T7 

1 d2n ’ lp (An) 2 

_ CO \ Pn d n np A p 

k ' m dn n m 

drinp A p 

2m dnl P > 

^ np 

2m dn n dn p ^ E ) 

dn n m 

1 (Aral 2 

T7 

i d2n "p(\„) 2 

dn np A p 

co i Pp dn np A p 
k ' m drip m 

2m dn n dn p ^ E ) 

11 PP 

2m dnl (/X P> 

drip m 


\ ( 5n n \ 
8n p 
Sp n 

V 5 PP ) 


= 0 . 


(36) 


The mode frequencies are given by the condition that the de¬ 
terminant of the matrix in Eq. ([36]) vanish. This result is a 


generalization to allow for entrainment of the results of An- 
dersson. Comer, & Prix ( |2004| i. Equation ( [36] ) illustrates that 
fact that, in nonlinear problems, density derivatives of n np oc¬ 
cur, as well as the quantity itself. 


5. APPLICATIONS TO THE OUTER CORE 
5.1. Equation of state 

The equation of state that we use is based on chiral effec¬ 
tive field theory (EFT), in which the symmetries associated 
with QCD are built into an effective Hamiltonian for nucleons 
(Epelbaum, Hammer & Meil3ner|2009 (.The parameters of the 



















































5 


theory are determined from nucleon-nucleon scattering and 
other low-energy nuclear data. The particular version of the 
theory that we shall use is that of Hebeler et al. (2013), in 
which an analytic fit is made to calculations for pure neutron 
matter and symmetric nuclear matter and an interpolation is 
made for proton fractions x = n p /n intermediate between the 
two proton fractions x = 0 and x =1/2 for which microscopic 
calculations have been made. Here 


n = n p + n„. (37) 

is the total density of nucleons. The nuclear part of the en¬ 


ergy per particle (without electrons) is given by Hebeler et al. 
( |2013l ) 


£ 

£o 


3 

5 L' 


r^/3 


+ (1- 


t)5/3 ] ® 


2/3 

— \ n 2 ! 3 


(38) 


-[aix(l-x) + a 2 ] — + y{nix{l-x) + r\ 2 ] — 


n o 


ho 


and the values of the parameters are Eq = 36.84 MeV, and 
cci = 6.14, rji = 4.02, a^_ = 1.4, and 772 = 0.9. 

This form is expected to be a reasonable approximation for 
baryon densities in the range ~ 0.08 — 0.2fm \ The energy 
density is the sum of the nucleon energy density and the elec¬ 
tron contribution 

E = ne+E e . (39) 


In the formalism described above, it is assumed that the num¬ 
ber of neutrons and the number of protons are conserved. This 
is a good approximation when the time scales of interest in 
the motions are short compared with the time scale for weak 
interactions. We have made no assumption about the ratio 
of neutrons to protons, but in the numerical calculations we 
shall concentrate on the case of matter in beta equilibrium, 
which should be a good first approximation for most of the 
life of a neutron star. The condition for beta equilibri um is 
that nipC 2 + dE/dn p + fl e = m n c 2 + dE/dn„ (Baym, Bethe, 
& Pethic k[l97lf which, with the neglect of the difference be- 
tween the neutron and proton masses, gives 


de/dx + fi e ss 0, 


(40) 


where fi e = dE e /dn e is the electron chemical potential, which 
for ultrarelativistic degenerate electrons is 

He = hc(37i 2 n e ) 1 / 3 . (41) 

Bulk matter is electrically neutral and therefore 

n e =n p . (42) 


The equilibrium value of the proton fraction calculated from 
Equations ( |39[ > and ( [40] ) is shown in Figure 1(a). For conve¬ 
nience, nucleon densities for matter in beta equilibrium are 
plotted as functions of baryon density in Figure 1(b). 



n [fm 3 ] 

Fig. 1.— Equilibrium proton fra ction and nucle on number densities calcu¬ 
lated from the equation of state of |Hebeler et al. 1(2013). (a) Proton fraction 
in beta equ ilibr ium as a function of nucleon number density, calculated from 
Equations {39} and {40} . (b) Equilibrium values of the nucleon densities n n 
and n p . Also shown are results for n np calculated from Landau Fermi-liquid 
theory with the Skyrme interaction SLy4 (see Eq. m 


5.2. Thermodynamic derivatives 
The second derivatives of the total energy density E, 


Eafi — 


d 2 E 

dn a dnp ’ 


(43) 


determine observable properties such as sound speeds. From 
Eq. ( [39] ) it follows that 



d 2 (tie) 

dn 2 

d 2 (we) 

dn 2 


Ep P — 


d 2 (we) 
dn p dn„ 


djle 

dn e ’ 


(44) 

(45) 

(46) 


We express derivatives with respect to particle density in 
terms of the variables n and x by using the relationships 


d 

dn„ 


_d_ 

dn„ 


+ and 
dx dn 

(47) 

1 — x d d 

n dx ^ dn 

(48) 


The results for the derivatives E a p obtained from Eqs. (44 »— 

S and (391 are plotted in Fig. 2. The quantity E pp has con- 
ttions from both protons and electrons, and we show the 
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Fig. 2.— Thermodynamic derivatives E a p and dfl e jdn e , the electronic 
contribution to Ep,,, for baryon densities in the outer core. The equation of 
state is taken from |Hebeler et aL](2013) . 

difference between E pp and the contribution from electrons, 
which in the absence of screening is d[i e jdn e . One sees that 
the electronic contribution to E pp is dominant. 



Fig. 3.— Sound sp eeds v in the absence of counterflow in units of c ob¬ 
tained from Equation (28) , as functions of the baryon number density (solid 
lines). The dotted lines correspond lines correspond to the velo citie s in the 
absence of coupling betw een neutrons and protons. Equation (30) (lower 
curve) and Equation ED (upper curve). The dot-dashed lines (v±) show 
the results in the absence of entrainment (det« a ^ = n„n p ). The modes corre¬ 
sponding to the three uppermost curves are dominated by motion of charged 
particles, while motion of neutrons is predominant in the modes correspond¬ 
ing to the three lowermost curves. 


5.3. Entrainment 

In addition to interactions between the densities of the vari¬ 
ous components, there are also interactions between the flows 
of the two components, which are reflected in non-zero val¬ 
ues of n np , an effect often referred to as entrainment. In the 
outer core of neutron stars, pairing gaps are expected to be of 
order 1 MeV or less, while nucleon Fermi energies are one 
or two orders of magnitude larger. Thus pairing contributes 
little to the total energy, and one may use Landau’s theory of 
normal Fermi liquids to calculate n np and Borumand, Joynt, 
|& Kluzniak| ( |1996] l find 


l np 


K Fp K F 

97Z 4 


-ft 1 


(49) 


where kp a are the Fermi wave numbers of neutrons and pro¬ 
tons, and f" p is the l = 1 component of the Landau parame¬ 
ter for the interaction between neutrons and protons. A gen¬ 
eral treatment of entrainment at nonzero temperature has been 
given by Gusakov & Haensel ( 2005) >. 

Most microscopic calculations of Landau parameters for 
nuclear matter have been performed for either symmetric nu¬ 
clear matter or for pure neutron matt er (For recent exam¬ 
ples see, e.g., Gambacurta, Lombardo, & Zuo| ((2011|; Holt, 
|Kaiser, & Weise| ( j2013| l),' and there is a need for further 
study of matter with proton fractions 5% of interest for 
neutron star cores. An exception is the work of Charnel & 


Haensel 
and ma< 


d|( |2006] >, 
de specif 


who gave a general treatment of entrainment 


specific calculations for effective interactions of the 
Skyrme type. For the standard form of the Skyrme interaction 
(Charnel & Haensel 2006j Eq. (23)), the entrainment comes 
solely from the terms involving gradients of the wave function 
and by direct calculation one finds 


/r=- 


k-Fnk-Fp 


1 + 2 Xl j ~ >rt2 f 1 2 X 2 


and therefore, from Eqs. ( |49[ and ( |50[ i, 

ttpp = 0 !nptlptlp 

in the notation of|Chamel & Haensel[(|2006]l, with 


m 

®np = 


tl ( 1 + —X\ j + ?2 f 1 + ~j*2 


(50) 


(51) 


(52) 


For the Skyrme interaction SLy4 developed especially for as- 
trophysical applications, t\ =486.82 MeV fm 5 , xi = — 0.344, 
t 2 = — 546.39 MeV fm 5 , and X 2 = —1.000 ( Chabanat et al.| 
|1998| > and therefore 


l np ' 


—1.567 fm i n„n„. 


(53) 


al. 

the 


As Charnel & Haensel (2006) showed, the 27 Skyrme interac¬ 
tions reco mmended for astrophysical applications by Stone et 
1998) give values for a np between 0 and —10.4 fm 5 , while 
Skyrme interactions developed by the Lyon group lead to 
values of around —1.5 fm 3 , with the exception of SLy230a, 
for which it is essentially zero. The wide range of values of 
n np that Skyrme interactions predict underscores the need to 
pin down its value better from more fundamental considera¬ 
tions. 

The conditions for stability to counterflow of the two com¬ 
ponents are that n, m , n pp , and det [n a p] are positive. If the third 
condition and one of the first two are satisfied, the remaining 
condition holds automatically. For the Skyrme interactions 
that have been considered above, n np is negative and therefore 
from Eqs. ( fl9| and ( [20] ) it follows that the first two conditions 
hold. Since 

det[n aj g] = ( n n n np )(n p n np ) n pp = n n n p n np n > 0, 

(54) 

the third condition also holds and matter is stable to counter¬ 
flow. 

5.4. Collective mode frequencies 

In Figure 3 we show results for the velocities of longitudinal 
collective modes. The velocities of modes in the absence of 
coupling between neutrons and protons are given by Eqs. 
and (|3T1( and the thermodynamic derivatives are taken from 
The dashed lines v± include the effects of E np but 
ects of entrainment are neglected (n, m = n„, n pp = n p , 
and n np = 0). Finally, the full lines include both the effects of 



nonzero E np and entrainment, Eq. (28 i. Entrainment affects 


the charged particle mode more than the neutron one since 
/n p is more than an order of magnitude larger than n np /n„. 


' L np 


The hybridization of the charged particle and neutron modes 
is relatively weak. When E np is nonzero but the effects of en¬ 
trainment are neglected, the velocity of the charged particle 
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Fig. 4.— Thermodynamic derivatives E a p in the inner crust and outer core. 
The equations of state (EOS) used are indicated in the figure. 



Fig. 5.— Speed of longitudinal sound-like modes in the inner crust and 
outer core of a neutron star as a function of baryon density. For the outer 
core, the results correspond to those in Figure 4. In the calculations for 
the inner crust, the neutron superfluid density was taken to be the density 
of neutrons outsi de nuclei. For the crust, the results are those of Kobyakov 
& Pethick 1 2016). The neutron mode (green line) tends to zero at the neutron 
drip density «,v/>« 2.2 x 10 4 fin \ The left vertical dotted line corresponds 
to the maximum density for which th e assum ption of spherical nuclei in the 
Lattimer-Swesty model is still valid ijKobyakov & Pethick|2013(. The right 
vertical line corresponds to the lower bound of density of unifo rm nuclear 
matter, below which it is unstable to fonnation of a density wave JHebeler et 
|al.|2013| . 

mode is raised, while that of the neutrons is lowered. Entrain¬ 
ment has little effect on the velocity of the neutron mode but 
further raises that of the charged-particle mode. 


It is instructive to compare properties of the outer core 
with those of the inner crust, where the protons reside in nu- 
clei. For the inner crust, values are taken from Kobyakov &| 


Pethick (2016 

i, which corrected a coding error in the paper of 

Kobyakov & 

Pethick (2013). These were based on the equa- 

tion of state of Lattimer & Swesty (1991). In Figure 4 we 


plot values of the thermodynamic derivatives E a p for the in¬ 
ner crust and the outer core. Despite the different equations of 
state in the cmst and the core regions, the values of E a p are 
rather similar at the crust-core boundary. 

Sound speeds across the crust-core transition region in the 
star are shown in Figure 5. It is interesting to note that, while 
the E a p are almost continuous between the crust and the core, 
the sound speeds exhibit significant jumps. At the boundary, 
the charged-particle mode is about three times slower in the 
crust than in the core; this is due to the fact that entrainment in 
the cmst is much greater than in the core by about one order of 
magnitude because in nuclei the number of neutrons entrained 
by a single proton is of order the neutron to proton ratio in 
nuclei, ~ 10, at the inner boundary of the cmst. 

6. DISCUSSION 


In this paper we have generalized to a two-component fluid 
Euler’s equation for a single component. The approach we 
have adopted is based on the Josephson equation for the 
phases of the condensate wave functions of the nucleons and 
the continuity equations. This makes possible a direct deriva¬ 
tion of the basic results. The nonlinear terms in the Euler-like 
equations have contributions proportional to density deriva¬ 
tives of the strength of the entrainment. These contributions 
do not affect small oscillations about a state in which the two 
fluids are at rest, but they do enter in, e.g., the condition for the 
two-stream instability. These terms are implicit in the work of 
) and they arise from the effects of entrainment 
chemical potentials]^] 

In some earlier treatments, the energy due to entrainment 
was regarded as part of an “internal energy” defined as the dif¬ 
ference between the total energy and the kinetic energy in the 
absence of entrainment ( |Prix|2004||Andersson, Comer, & Prix] 
|2004| , but the approach presented here shows that it is natural 
to treat the energy due to entrainment as part of the “kinetic 
energy”. In this way it is made clear that the nucleon chemi¬ 
cal potentials contain contributions proportional to derivatives 
of the entrainment energy density with respect to the neutron 
and proton densities. The thermodynamic potential appropri¬ 
ate when the system is specified by the number densities of 
neutrons and protons and the phases of the condensates is the 
Hamiltonian, that is the total energy of the system, while its 
Legendre transform, 

S=p„-j„+p P -jp-£' tot 

= ^X/»( n_1 )aj8 ia-jp-E, (55) 

Z ap 


Mendell (1991 


on the nucleon 


is the potential appropriate when the current densities are re¬ 
garded as the variables. Here the matrix n 1 is the inverse of 
the matrix with elements n a p. Numerically, the first term on 
the right side of Eq. ( |55j ) is equal to the kinetic energy, Eq. d9]). 

Our calculations sKow that in the generalizations of Eu¬ 
ler’s equation to two-component superfluid hydrodynamics 
first and second density derivatives of the entrainment func¬ 
tion n np appear. Nonlinear effects in superfluid hydrodynam¬ 
ics have been investigated in a number of different contexts 


2004; Gusakov & Andersson 

2006; 

Glampedakis, Andersson, 

& Samuels son 

2011, Haske 

1 201 

; Link [2*0l2f|Passamonti 


& Lander|2012 


, and an important task for future work is to 
investigate to what extent results are altered by the nonlin¬ 
ear terms derived in the present article. It is also necessary 
to reexamine how the terms obtained from a Hamiltonian ap¬ 
proach are reflected in the Lagrangian and hybrid approaches 
used in other work. 

In this article, we have assumed that the flow is irrotational, 
in the sense that V x p„ and V X p p vanish. We leave for 
future work the incorporation of electromagnetic fields, vor¬ 
tices, and rotating frames of reference. An additional direc¬ 
tion for investigation is the effect of nonzero temperature, 
which results in the appearance of a normal fluid of excita¬ 
tions. 

As applications, we have considered oscillations of uni¬ 
form neutron star matter. We have generalized the treatment 


9 The original work of Andreev & Bashkin (1975} on entrainment in the 
helium liquids did not mention explicitly the entrainment contributions to 
the chemical potentials. However, this had no influence on the applications 
described in that paper, which were to linear modes. 
























































8 


of the two-stream instability given by Andersson, Comer, & 
|Prix| ( |20(j4] l. To make realistic estimates of the conditions un¬ 
der which the two-stream instability can occur in neutron star 
cores, it is necessary to take into account damping: in particu¬ 
lar, it is important to include pair-breaking processes that will 
set in at wave numbers of approximately A /vp ~ kp(A/Ep), 
where A is the superfluid gap of a component, vp its Fermi 
velocity, kp its Fermi wave number, and Ep its Fermi energy. 
These wave numbers are much less than the respective Fermi 
wave numbers. 

Velocities of sound-like modes in the outer core in the ab¬ 
sence of counterflow have been calculated. In particular, we 
have generalized the discussion of|Bedaque & Reddy|(|2014 1 
to allow for entrainment, and we have used recent calculations 
of the equation of state to evaluate the thermodynamic deriva¬ 
tives. Extensions of this work to shorter wavelengths and to 


calculate damping of modes by the electrons will be reported 
elsewhere pCobyakov et al.|2016|l. 
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APPENDIX 

COMPARISON WITH EARLIER WORK 


The results of the present work agree with t he work of Mendell ( 1991). He re we compare these results with those of other 
studies that are also based on Mendell’s work. In jAndersson & Co mer (2001) and Anders son, Comer, & Prix| ( |2004[ > the Euler-like 
equation for the neutrons has the form 


dt dx. 


r , 1 d{l^ CP <3v„; 

[l ’ni T £/( \}’pi I'm)] 5 ^ ~ f £n\}’pj Di j) -j y — 0, 


where the indices i and j refer to Cartesian coordinates, and the index j = 1.2,3 is to he summed over. Here 

-l 


£« — 


2a 


ftpYlnp 


L np 


det n 


a/3 


1 - 


L np 

Tl n Tlp 


(Al) 


(A2) 


and we denote by p^ CP the chemical potential used in those papers, which is different from those employed in the present article. 
Comparis on w ith the present work is simplified by observing that the combination v„ + e„(\ p — v„) is what we denote by p ,,/m. 
Equation ( |A1[ ) may therefore be written as 

Pnj d \ p ni d fl£ CP d d p m 

mdxj)m + dxi m + £n{Vpj V ' ,j) dx^ En{Vpj Vnj) " (A3) 


d_ 

dt 


dxj m 


If the problem is one-dimensional, with all variations in the x-direction, one finds 


d_ 

dt 


Pnx r > 

m dx 


Pnx 

m 


d_ 

dx 


I, ACP i 

^-2**-v„n=o. 


(A4) 
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We now contrast this result with the one found from the present work. Equation (22 1 for the one-dimensional case reads 

2 dn„ ^ 


( d 

! Prvc <5 ^ 

Pnx 

+ T X 

( Pn 

i 

n„n p 

\dt 

m d xj 

m 

l m 

2 

detrc a p 


L np 


iy px ^nx) 1 — o, 


where we have made use of the relation 


decaff p p -p„ 

n n n p m 


(A5) 


(A6) 


which follows from Eqs. 0-0 and ( [TT] ). The terms containing v„ — v ; , in equations ( |A4[ > and ( |A5| do not agree. We have been 
unable to find in the literature an explicit expression for p^ cp . If it is to be identified with the chemical potential in the absence of 
flows (what we denote by p n ), there is a conflict. There is too if it is identified with p„ — (dn n p/dn n )(p p — p n ) 2 /2m, the chemical 
potential with the entrainment contribution but not that from the flow of the neutrons. Similar conclusions apply for the protons. 









